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Abstract 

Although aircraft have been routinely flown at transonic speeds 
for the last two decades, the designers of these machines have had 
to formulate their designs almost exclusively on the basis of ex- 
perience, In contrast to subsonic or supersonic aircraft design in 
which a welter of analytical, experimental, and numerical techniques 
exist. Since this type of fluid flow is characterized by complex 
viscld-inviscld Interactions, the development of fast numerical models 
of the full Navler-Stokes equations has promised to alleviate this 
situation. One such model Is the approximate factorization algorithm 
Introduced by beam and Warming and implemented by Steger, et al. 

This research, then, is principally concerned with an Independent im- 
plementation of this numerical algorithm and initial studies of its 
ability to efficiently and accurately describe transonic flow about an 
NACA 64A010 airfoil section. 

The approximate factorization algorithm is developed from the non- 
dimensional, conservative, vectorized Navier-Stokes equations expressed 
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In curvilinear coordinates. Equations of state and transport coef- 
ficient relations appropriate to atmospheric air are appended to close 
the system of partial differential equations. An algebraic turbulence 
model due to Baldwin end Lomsx is also incorporated Into the equation 
set. The coordinate generation mefl'od developed by Thompson, et al is 
used to produce the desired coordinate transformations. Boundary con- 
ditions on the airfoil surface are formulated so as to allow suction 
and/or blowing from the surface and to emulate either an isothermal 
or adiabatic wall. Outer boundaries are placed ten chord lengt’is from 
the airfoil and their boundary conditions formulated so that the in- 
flow properties can be varied and the outflow properties determined 
by extrapolation. Fourth-order artificial viscosity proves to be nec- 
essary for high Reynolds number flows. 

This algorithm was verified by investigating the flow about an 
NACA 6AA010 airfoil at 0“, 2®, and 3.5° angle of attack for free-stream 
conditions of 2 x 10^ Reynolds number and 0.8 Mach number. The flow 
was Initiated by either gradually decreasing the degree of fluid pen- 
etration of the airfoil from total to none, or by using a body force 
to gradually accelerate the airfoil and its attached coordinate system 
from zero velocity to frce-stream values. One zero degree case was 
run as laminar flow while all the other cases used the algebraic turbu- 
lence model. Also, cases were run to verify the model's ability to main- 
tain a free-stream solution. As an aid in evaluating the results, a 
set of Mach number contour plots and coefficient of pressure graphs 
were prepared. 
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Overall reaulta were In good qualitative agreement with the wind 
tunnel data aeta of Johnaon and Bachalo. Unfortunately, while non- 
dlmcnalonal times of alx were attained, numerical difficulties pre- 
vented any case from reaching a true steady state. In the last test 
case attempted, that at an angle of attack of 3.5°, there was no doubt 
that a shock was forming on the airfoil ana that separation had occurred. 
Computer times for the 113 x 51 grid used were encouraging, averaging 
35 eeconds/tlme-Btep on a Unlvac 1100/80 computer and 5 seconds/time- 
step in scalar mode on the Cyber 203 (Star) . It Is concluded that ap- 
proximate factorization techniques, while they still need some work, 
can definitely be used to advantage in at least two-dimensional tran- 
sonic flow problems. 
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I. INTRODUCTION 


Although the art of modeling fluid flows has progressed to the 
point where important parameters of certain complex flows can be 
routinely and confidently predicted, there are still some inter- 
esting flows for which current methods of analysis (experimental, 
analytical, or numerical) give approximate results at best. One of 
these difficult fluid flow problems is the description of transonic 
flow over geometries of practical interest. In particular single- and 
multi-element airfoils. This type of flow is of particular interest 
to the numerical analyst since it is of major concern to the aircraft 
industry and since it proves to be inherently too complicated a flow 
to be fully described by either pure analysis or experimental tech- 
niques. Until recently, the same could be said of the available nu- 
merical methods, however, with the appearance of improved ADI (alter- 
nating direction, implicit) techniques, this is no longer the case. 

The major thrust of the research reported herein was to establish Just 
this premise and to verify the numerical algorithm developed for this 
purpose . 

The numerical algorithm Itself is based upon a finite difference 
representation of the full, Reynolds-averaged, Navier-Stokes equations 
as developed by Beam and Warming [c] and including some of the implemen- 
tation ideas due to Steger [b]. The set of equations is closed by 

using the standard approximations for atmospheric air and including an 
algebraic turbulence model developed by Balwin and Lomax [a] from the 

work of Cebeci [p]. Chief among the attributes of this algorithm is 

that it is a non-iterative, second-order accurate, implicit formulation 


of th« conservative form of the governing equations which remains stable 
and accutate for Courant numbers much larger than unity and which at- 
tains a steady-state solution Indeopadcat of the time step. This 
method offers the practical advantage of being able to calculate ac- 
curate transient and steady-state solutions In time periods measured In 
minutes, rather than the hours of older methods, at the cost of re- 
quiring large amounts of data storage even for two-dimensional problems. 

Almost all recent numerical models of fluid flows are based upon 
forms of the Navler-Stokes equations expressed In curvilinear coordi- 
nates. That this should be so Is apparent from the fact that most 
fluid flows of practical Interest have Inherently disparate length 
scales arising from the geometry and/or the equations themselves. Thus, 
curvilinear coordinates offer an automatic adjustment to these varying 
length scales and the bonus of simple boundary condition specification. 
One of the most widely used techniques for automatically generating 
such coordinates, that due to Thompson, Thames, and Mastln [d], has been 
adapted to this research. It allows almost complete control over boun- 
dary location and provides for concentration of coordinate lines about 
the airfoil. 

Of course any new numerical algorithm and in particular any com- 
puter coded version of any numerical algorithm, needs to be extensively 
verified against known solutions and/or trusted data sets. Although, 
as mentioned in the opening paragraph, there are scant theoretical so- 
lutions for the flows considered here, there have been some recent ef- 
forts at developing a systematic set of experimental data for use in 
testing numerical models. The paper of Johnson and Bachalo [m], detail- 
ing wind tunnel test results for the NACA 64A010 airfoil section, was 
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p^irtlcularly influential In choosing the flow c editions to be modeled. 
Their results Include Inter ferograms, calculated Mach number contours, 
and an extensive compilation of turbulence parameters for a freestream 
Mach number of 0.8, and Reynolds number of two million, for various 
angles of attack. 

The dissertation contained In this volume basically follows the 
same logical de^ ? ..opment as the foregoing discussion. Following a sur- 
vey of the literature which this author considers to be germane to the 
topics covered, the nondlmenslonal, Reynolds-averaged, vectorized, 
fully conservative Navier-Stokes equations are put Into a form appro- 
priate to curvilinear coordinates. Also, the equations of state and 
the functional form of the transport coefficients are presented as part 
of the development. The succeeding chapter details the transformation 
of these partial differential equations into a set of difference 
equations which form an ADI scheme Involving a non-iterative, block- 
tridiagonal inversion during each sweep. Determination of the trans- 
formed coordinates and appropriate finite difference boundary conditions 
are also included In this chapter. Next, the verification trials for 
the computer coded algorithm are set forth. Results for the NACA 
6AA010 airfoil at 0°, 2®, and 3.5° angle of attack under transonic con- 
ditions are presented and checked for internal consistency. The avail- 
able experimental data is also compared with these results and appro- 
priate comments made about their degree of and/or lack of correspondence. 
Lastly, in the concluding chapter, some general and some specific con- 
clusions about the algorithm and data set are made and some recommenda- 
tions for future research work given. 
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It If. hoped that the casual reader, should 
will find the exposition clear and enlightening 
student will find the presentation to be useful 
stimulating. 


there be such a person, 
and that the serious 
self-contained, and 
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II. LITERATURE SURVEY 


The purpose of this chapter is to acknowledge those authors whose 
published work has contributed to this research in some way; it is also 
Intended to provide the Interested reader with some related references 
outside of the narrow scope of this research. However, in no way is 
this chapter to be considered as a definitive, exhaustive or complete 
treatment of the significant publications dealing with the topics 
touched upon herein. 

Although some current analytical work is being done on the problem 
of transonic flow, most of it is based upon the techniques developed 
earlier in this century before the advent of electronic computation. 
Landau and Llfshitz [q] in their survey of the state of fluid mechanics 
(1959) give a rigorous introduction to the analysis of transonic flows. 

Some of the earliest numerical work in this area was based on 
various simplifications of the full potential flow equations; indeed, 
this approach appears to be a well-entrenched field of study for at 
least the foreseeable future. A good source book for some of the cur- 
rent (1978) potential methods is listed as reference [r]. 

Numerical solutions of the full Navler-Stokes equations became 
practical for transonic flows when MacCormack [s] introduced his 
"rapid solver" algorithm. However, being an explicit technique, it po- 
ssesses certain relaxation limitations. Although implic'w algorithms 
have a long history of development in this field, it was the appear- 
ance of approximate factorization type techniques which made them more 
than of just theoretical interest. These algorithms were first advanced 
by Beam and Warming [c,f] and by Briley and McDonald [t] independently. 
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The B««m and Warming formulation has been successfully applied to a 
variety of transonic flows by Steger [b], Pulliam and Steger [g], and 
Steger and Bailey [n]. 

The advantages Inherent in the use of curvilinear coordinates 
have been made readily available by the pioneering work of Thompson, 
Thames, and Mastln [d]. Recent modifications of the original technique 
are detailed in reference [e,J,n]. 

Experimental data sets which are particularly adapted to verifica- 
tion of transonic algorithms have been developed in recent times. In 
addition to the NACA 64A010 airfoil data used in this research (Johnson 
and Bachalo, reference [m]), the data set gathered by Seegmlller, 
Marvin, and Levy [u,v] for an 18% thick circular arc airfoil are to be 
reconmended . 
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1X1. VECTORIZED NAVIER-STOKES EQUATIONS IN CURVILINEAR CXX)RDINATES 


Although many almpllflad versions of the Nsvler-Stokes equations 
have been successfully applied to various transonic problems, In 
general this very complex type of fluid flow can only be adequately 
described by the full set of equations. However, due to certain un- 
solved problems In turbulence modeling, constraints Imposed by the 
numerical technique used, and the desirability of keeping the analysis 
as simple as possible, a number of assumptions and restrictions are 
required. For example: due to the Impossibility or Impracticability 

of using enough grid lines to resolve the small scale eddies In turbu- 
lent flow, the Reynolds averaged Navler-Stokes equations are used along 
with an algebraic eddy viscosity turbulence model. This chapter, then, 
will present the particular form of the Navler-Stokes equations, the 
equations of state, the constitutive equations, and other subsidiary 
relations which were chosen to model the transonic flow problems dis- 
cussed In this paper. 

A. Problem Description In Cartesian Coordinates 

The governing differential equations used to model the transonic 
flow of air In this work are the two-dimensional, Reynolds averaged 
Navler-Stokes equations for a Newtonian fluid which obeys a Fourier 
heat conduction law: This system of four partial differential equations 

(continuity, x-momentum, y-raomentura, and energy equations) can be ex- 
pressed in the following vectorized, nondimenslonal form: 
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where 



'p ■ 
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■ m 

pv 
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puv 
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E(q.q^tqy) 


■ 

0 

yVx + yUy 

Xu + (X+2y)v 
X y 

yuv + Xvu + (X-f2y)w + yuu + -;r" T 
X X y yPry 


(3.4) 


In addition, since air can be assumed to be a (thermally and calorically) 
perfect gas for the flow studies, the following nondlmenslonal equations 
of state were used: 


T “ yC^ - |(u2 + y2)] 

P - (Y-l)[e - j p(u^ + v2)] 


(3.5) 

(3.6) 


The transport coefficients were obtained by assuming that the Prandtl 
number is effectively constant, that a Sutherland viscosity law is 
valid, and that Stokes law can be used. These relations are listed 
below in nondlmenslonal form: 
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(3.7) 

(3.8) 


k ■ )j 

W - [1 + (y-1)m2Sj] (Y-l)Sl 
X - - I M (3.9) 

Th« no- idlnena tonal quantltlaa appaaring in equations 3. 2-3. 9 are de- 
fined in the List of Symbols following the Table of Contents; details 
of the method of nondimensionalization are contained in Appendi)' A. 

An algebraic turbulence model can be incorporated into the above 
set of equations by multiplying the viscosity coefficient, v , by the 
factor: 

1 + (3.10) 


except where it replaces the thermal conductivity, k , (equation 3.7); 
then the factor 

1 -f ~ ^ (3.11) 

Pr^ p 


must be used. The eddy viscosity, , and the turbulent Prandtl 
number, Pr^ , depend upon the turbulence model used. Due to Its rel- 
ative simplicity and Its prior use in transonic flow calculations, the 
Baldwln-Lomax [a] algebraic turbulence model was used when required. 

Of course, no problem is fully specified until the initial and 
boundary conditions are stated. However, since this research involved 
distinct problem types (l.e. couette flow and transonic flow over 
multi-element airfoils) and since the particular boundary and/or initial 
conditions used are rather intimately related to the numerical solution 
procedure, further discussion of this topic will be deferred. 
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B. Problam Dtscriptlon in Curvllinaar Coordlruf 


Th« conq>utaclonal grid on which this act of partial dlfferantial 
aquatlona ara aolved uaually doaa not fom a Cartaalan coordlnata 
ayataa. Thua It la advantagaoua to ra-axpraaa the problem In tarma of 
more general curvilinear coordlnatea while retaining the atrong con- 
aervatlon form of the Navlar-Stokea aquatlona. In the notation com* 
monly choaen (aea Stager [b], for example) the following coordinate 


transformation la defined: 


T ■ t 

(3.12) 

i - C(x,y,t) 

(3.13) 

n ■ n(x,y,t) 

(3.14) 


Using this transformation, equations 3. 1*3.4 can be written as: 


3t 


+ A*(q*) + h 


i. 

3n 


, J_rl_ 

Re**3C 






where 



(3.15) 

_* _ 1 „ 
q - j q 

(3.16) 

A’'(q*) - jCt^q + 

(3.17) 

B*(q*) - j[nj.q + n^A+riyB] 

(3.18) 


(3.19) 

E*(q*.<i*5 .q*> ■ + 

(3.20) 


The important details of this transformation are contained in Appendix 
B. In that, for the most part, the remainder of this work will deal 
with this generalized vector equation, the superscript will be 
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■upprassdd axctpt wh«r« needed for clarity. 

The equaclone of etate (3.5 end 3.6) end the trenaport coefficiente 
(3. 7,3. 8, end 3.9) ere not effected by a coordinate trena format ion, 
except that the dependent verieblea muat now be interpreted ea functione 
of C,n, end t (i.e. P(C,n,T) = P(x(t,n,T) ,y(C*n,x) ,t(x))) . 

Since the numerical algorithm which wee eelected to eolve equation 
3.15 muat treat croaa-derlvetivea apecielly, it ie necessary to split 
the viscous vectors D and E into vectors which contain only K- or only 
0 - derivatives. Thus equation 3.15 is written as: 


It + Ic *<■') ♦ It "«» ■ 


where 



m s 

P 


pU 


■pv 

1 

pu 

pv 


puU + t P 

X 

pvU + C P 

y 

’ "-I 

puV + n P 

X 

: vV + n P 

y 


e 

to m 


(e+P)U-CjP 


(e+P)V-n^P 


_1 

J 


°2 “ J 


®1^ ^ Ve 
Vc + 


a^uu^ + a2<vu^ + uv^) + 

m ■■ 

0 


b,u + b-v 
In 2 n 

b^uu^ + b3vu^ + b2uv^ + b^w^ + b3l^ 


(3.21) 


(3.22) 


(3.23) 


(3.24) 
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0 


"l- J 


b,uu^ + + b3uv^ + b^w^ + bjT^ 


(3.25) 


I 

J 


c.u + c.v 
In 2 n 

c_u + c.v 
2 n 3 n 

c.uu + c,(vu + uv ) + c.vv + C.T 
In 2 n n 3 n n 


(3.26) 


In equations 3.15-3.26, the transformation Jacobian, J , is defined 


as: 


J - C,ny - Cyn^ (3.27) 

The contravariant velocities, U and V, are given by: 

u • Ct + V 

V - nt + n^jU + TiyV (3.29) 

Also, the viscous coefficients are defined as: 

a - (X+2y)?2 + 

1 X y 

32 - (X+U)^^5y 

®3 " 

“4 ■ 


12 


b 

b 

b 

b 

b 


1 

2 

3 

4 

5 


*^1 


(A+2li)Cj^njJ + UCyHy 

X5 n + Mt n 

X y y X 

u5 n + AC n 
X y y X 

wC n + (A+2y)c n 

X X y y 

Pr^^x\ ■*■ 


(A+2p)n5 + mo5 

X y 

(A+u)n n 
X y 


wn,^ + (A+2p)n2 

■^(n^ + n^) 

Pr X y 


(3.30b) 


(3.30c) 


Note that except for the doubling of the number of viscous vectors, 
the curvilinear equations (3.21-3.26) are not much more complex than 
the Cartesian equations (3. 1-3.4). This can be very useful in that 
some relations developed in Cartesian coordinates may be directly con- 
verted to curvilinear coordinates by maintaining the proper correspond- 
ence of terms. 

Since none of the problems Investigated in this research involved 
the use of a time-dependent coordinate system, henceforth, the deriv- 
atives and will be taken to be zero. 

As was the case with the equations of state and the transport co- 
efficients equations, any boundary or initial conditions which do not 
involve spatial derivatives will remain Invariant in form. However, 
there are two common types of boundary conditions which are derivative 
relations. One of these is the "adiabatic wall" assumption, or more 


g«n«rally, a apaclfleatlon of tha haat flux, q" , through soma portion 
of tha boundary of the flow domain. In nondimanaional form this con> 
dition can ba expressed as 

which follows from tha Fourier haat conduction law and tha nondimens ion- 
alization (see Appendix A). Tha symbol "7^" indicates a derivative 
normal to tha boundary surface. In curvilinear coordinates this re- 
lation becomes: 

q" - ^ (oT^ + eT^)/^ (3.32a) 

where 

o«n^ + n^ , (3.32b) 

X y X X y y 

The other derivative boundary condition is one which arises from 
the fact that whereas the analytical formulation of a problem only 
allows three of the four dependent variables to be specified at a 
solid boundary, the numerical formulation requires all four to be spec- 
ified. To maintain some degree of consistency with the analytical prob- 
lem, this ’’extra" boundary condition should take the form of an extrapo- 
lation from the Interior of the flow field. Thus, one of the most 
connionly used additional boundary conditions is to specify the normal 
pressure derivative in accordance with either boundary layer theory or 
the momentum equations. Using the second approach, the momentum 
equations and continuity equation may bp combined c-o form: 
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V + t‘Vt + Vt’ * + Vi’" * ■ 

^ ([-.,(., 'uj + + n,(.^«5 + .Jv^)^] + 

[n^CbJu^ + bjv^)^ + „^(bju^ + biv^)j] + 

[.f,(bju5 + bjvj)^ + HyCbj'Oj + b-v^)^] + 

f"x«=l'“n * '2%>n * "y<'2“n * '3%’n” 


(3.33) 


where 

VjjP - (aP^ + 0P^)/»^ (3.34) 

and the primed coefficients are the viscous coefficients (3.30) 
divided by the Jacobian "J". Further details on the development of re- 
lation 3.33 are contained in Appendix C. 


C. Curvilinear Coordinate System Generation 

The use of a curvilinear coordinate system is desirable because, 
among other reasons, it allows easy application of boundary conditions 
since each boundary can be made to coincide with a coordinate line; 
it can concentrate coordinate lines in one region, place a minimum 
number in another region, and smoothly transition from one to the other; 
and it allows numerical algoritnms to be relatively independent of the 
particular geometry of a problem. However, for curvilinear coordinates 
to be truly useful, one must utilize a method of coordinate generation 
which, in addition to providing the advantages of the previous sentence, 
must be relatively simple to implement, must generate "smooth" (l.e. 
second order and higher order derivatives small In some sense) , single 
valued coordinates, and, at least for some of the applications 
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eonaldarcd In this paper, must be capable of handling multiple bodies 
easily. The method chosen for this work Is the ''boundary-fitted 


coordinate" algorithm due to Thompson, Thames, and Mastln [d,e], 
which achieves Its purpose by solving a Poisson equation In the curv- 
ilinear coordinates C and n> 

Briefly, this technique requires that the coordlnatea satisfy the 
elliptic system: 

“c 

5 +5 . . p (5 (3.35) 

XX 'yy j2 c 


n + 

XX 


yy 



Qg(f.,n) 


(3.36) 


where normally the boundary conditions to this system are specified so 

that the body(ies) lie along certain of the ri coordinate lines with 

some desired distribution of C lines terminating on them. The functions 

P and Q are specified so that the coordinate spacing in the Interior 
c c 

of the field is an approximation of the one desired. In order to ob- 
tain the calculation advantages of solving these equations on a curv- 
ilinear coordinate system, they are transformed so that it is the x and 
y distribution that is determined. That is, the following elliptic 
system is actually solved: 


a x__ - 2B x_ + Y X ■ -(a x P + Y_* Q„) (3.37) 

c c Cn c nn c c c n c 


a y, + Y y “ ”(» Ve-P + Y y Q ) (3.38) 

c^55 c'5n ’c'^nn c'C c 'c-'n c 
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In practice, the above ayatam la dlacretiaed and aolved by atandard 

methoda of numerical analyala. Further detalla on the development and 

uae of thla method of coordinate generation and the definition of the 

parametera a , 3 , y » and J are contained In Appendix D. 
c c c c 


IV. NUMERICAL FORMULATION 


It la an unfortunate fact that mathamatleal analyala of nonlinear 
pnrtlal differential aquatlona» euch aa the Navlar-Stokaa aquatlona» 
can not currently provide anything Ilka a cloaad form aolutlon to thaae 
equatlona for general boundary condltlona. Thue, though Itaelf more 
of an art than a science » methode of numerical analyala must of neces- 
sity be applied to these problems If more than a qualitative descrip- 
tion of the flow la desired. One of the more common numerical methods 
used to solve the Navler-Stokes equations la to approximate this set 
of partial differential equations by an equivalent set of difference 
equatlona In such a way that they are consistent with these differential 
equations. In addition, their solution algorithm must produce a se- 
quence of Intermediate Iterative values which converge to the actual 
solution to some estlmatable degree of accuracy. From a practical 
standpoint, the iterative algorithm should possess the property of a 
rapid rate of convergence, both numerically (small number of iterations) 
and computationally (small number of operations). Thus, one of the most 
important aspects of this research was the selection and development 
of the numerical algorithm used to solve the Navler-Stokes equations 
presented in the previous chapter. 

The difference scheme chosen for this work was the approximate 
factorization algorithm due to Beam and Warming [c,f]. In the form 
used in this research, this algorithm is an implicit, second-order ac- 
curate in time, non-iterative (in the sense that each time-step is cal- 
culated once only), unconditionally stable (in the linear approxi- 
mation), three-time-level scheme. Its "delta" formulation Insures 
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« 


« 


that I although tha croaa-darlvativa tarna ara traatad axplleltly, 
tha sehaina ramaina aacond-ordar accurata and unconditionally atabla; 
raquiraa only two-tlma-lavala of atoraga; and produeaa a ataady-atata 
solution which la Indapandant of tha tlma-atap alsat Tha noat attrac- 
tlva faatura of tha algorithm la that» due to Ita apatlal factorl- 
aatlon« It forms an ADI type of schema In which each sweep Involves 
the Inversion of a block-trldlagonal matrix. In principle this ap- 
proximate factorization algorithm promises to solve complex, two- 
dimensional fluid flow problems on current high-speed scientific com- 
puters within an economically reasonable time (l.e., ten to thirty 
minutes) . 

The remainder of this chapt will attempt to put the assertions 
of the previous paragraph into a more concrete form. First the delta 
formulation is presented along with an indication of the method of 
spatial differencing. Then the techniques employed to impose the boun- 
dary conditions are detailed. Lastly a brief discussion of the problem 
of forming the geometric coefficients is included. 


A. Beam and Warming Approximate Factorization Algorithm 

The most general, consistent, three-tlme-level, linear expression 
relating the conservation variable, q , and Its time-derivative, q^ , 
is [f]: 




. V • 




, n-1 


( 4 . 1 ) 


with truncation error of: 

+ 0-|)4t2q"^ + 0(At3) 


( 4 . 2 ) 
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t 

\ 


t'. 

I 



"A" It tht ututl forward difftrtnct optrttor, and 6,^,4* art 
arbitrary raal conatanta. Nota that tinea t • t> thay will ba utad 
intarehangaably in thia paper. Haneaforth, the paramatar ^ will alwaya 
ba taken to ba zero and* to naintain aaeond-ordar accuracy, the fol- 
lowing relation ia aatabliahad between ^ and 6: 

’■♦*5 „.3: 

Now the time-darivativaa of the conaervation vector are evaluated 
from equation (3.21) and aubatltuted into equation (4.1) with the 
reault : 


n A - a AD.+AD, „ . AE.+AE, „ 


D,+D 


E.+E, 


. AtrB, . . r‘^2.n . 8, „ . ‘'r‘'2xn., . 6 . n-1 

* Wh*-* * * “rT-> 3 + A'"' 


This equation is not suitable for direct use since the time-differenced 
vectors on the right-hand-side involve nonlinear functions of q”^^. 
However, due to special properties of the flux vectors, (A,B,D, and E) , 
under a suitable assumption on the transport coefficients, these dif- 
ferenced vectors can be expressed in terms of their Jacobian matrices 


as follows: 

Aa" - ~ a” Aq" - p"Aq" (4.5a) 

Ab" - " Q%“ (^.5b) 

4d" - [g|-Dj Aq"]^ - (WV)^ (4.5c) 

n 
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(4.5t) 


AeJ - [^eJ Aq"]^ - (Z%“)^ (4.5f) 

where Che croee-derlveClve netricee (AD 2 end AE^) heve been lagged In 
Cine In order Co feciliCeCe Che fecCoricaClon deecrlbed leCer. Theee 
expreaelona are exacC Co Che order of Ac^, and chue chair uee In 
equeClon (4.4) will noC degrade Che accuracy of che algorlchm. AlChough 
noC aCrlcCly neceaaery Co che developmenc, IC le convenlenC Co also use 
Che following exacC relaClonshlps becween che flux vecCora and Chelr 
Jacobian maCrlces: 

Dj - Wq^ , Dj - Xq^ 

h " » h " 

and 


Wq - Xq - Yq - Zq - 0 (4.7) 

A deCalled derlvaclon of Cheae equaclons (4.5, 4.6, and 4.7) la con- 
Calned In Appendix E; the expllclC form of the Jacobian matrices 
(P, Q, W, X, Y, and Z) are contained in Appendix F. 

Using relations (4.3) and (4.6), equation (4.4) may be re-expressed 


as: 


. 0Atr. 3 - I 3^ ..V . / 3 « 1 3^ » vi^n^.n 

“ * * ‘TS** • 2 > J> '“1 


A + ^(Wq, + Xq.)] + B + i(Vq, + Zq.)])" 


1+4) 3 c 


3'0‘ 


Re' 




(4.8) 


where "I" Is the four-by-four Identity matrix and the bracketed term on 
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th« left of Aq” !• undorstood to bo o dorivotlvo oporotor. Slnco 
oquotlon (4.8) io of order 4t^ end elnce ie itself of order At, 
this equation een be "epproxloetely factored" into the following sot 
of equations with no lose of accuracy: 




„ . 6 Atr a ^ 1 92 n TH 


^n+1 _ „n . .^n 
q ■ q + Aq 


where Aq“ is en intermediate velue and RHS (4.8) represents the rlght- 
hend-slde of equation (4.8). Although equation eat (4.9) repreeente 
the basic approximate factorisation algorithm, the spatial derivatives 
must also be discretised for it to be of practical utility. Using the 
central difference approximations developed in Appendix G, equation 
(4.9) can be put in the finite difference form: 

“*®1,J " ^^*1+1,J^‘*1+1,J+1 ■ ‘‘l+l.j-1^ " 


+ 2 [(R%" ^i+l,j+l " ^ 1 +l.J-l 




n._n-l, 


^i*ui ■ ““J.J 

(4.10b) 

S.J-I + ♦ «-«S,j+i <.i*i • 

(4.10c) 


"i.d • "i.j * '"i.j 


(4.10d) 


where R ■ -|(XfY) and the coef flclente ^ end ^ have been eb- 

aorbed by the netricee U, X» Y, Z end A* B» P, Q, reapectlvely. 

Thie final net of difference equatlonn (4.10) elao deflnee the 
bealc eolation algorithm. Flrat the explicit vector RHS la calculated 
over the Interior of the field. Then each n-llne la ewept and a 
block-trldlagonal loatrlx Inverted to determine the Intermediate valuee 
4q" .. Similarly* each C~Hne la ewept and the reaultlng blocfc-trldlag- 

*t j 

onal matrix Inverted* thus determining the difference conservation 
vector 4q" . Lastly the conservation vector q" . is determined from 

*»j *»j 

equation (4.10d). In essence* then, this method Is an ADI type of Im- 
plicit algorithm In which just one Iteration per time-step Is performed 
due to the direct Inversion of the block-trldlagonal matrices. For 
completeness* the block-trldlagonal Inversion algorithm is Included In 
Appendix H (see also Steger [b], his Appendix 1). 


B. Initial Conditions 

Even though it is the primary objective of most fluid flow analyses, 
Including this one* to achieve a solution which is independent of the 
Initial state of the fluid* thus making the boundary conditions all im- 
portant* the eminent difficulties Involved in getting a numerical sim- 
ulation started at high Reynolds and Mach numbers results in the initial 
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\ 

I 

i 

j 

conditiona taking on a eoaputatlonal Inportanca which can not ba | 

naglactad. Xndaad, tha atarting problan oftan dictataa* at laaat | 

aarly on» tha boundary conditiona to bo uaad and oftan nodifiao tha I 

i 

baaic algoritha Itaolf* In tha raaaarch raportad in thia voluna, thraa 
difforant aota of initial conditiona and thair aeconpanyfug modifi- 

cationa to tha boundary conditiona and algoritlun warn triad. | 

Tha oimplaat initial condition, at laaat from tha atandpoint of 
initial formulation, in tha ao-callad impulaiva atart. In thia caaa j 

I 

tha initial propartiaa of tha flow fiald arc tha aama aa tha fraa-atraam 
propertiao except for tha valocitiaa on any aolid boundariaa. Thia 
initial condition apaciflcation haa two major diaadvantagea, both due 

I 

to the uaually large velocity Jump near a body. Moat deleterioua ia ] 

• 1 

the fact that the algorithm muat be altered ao aa to rapidly diffuae | 

the Btrong velocity gradient at the body. One auch alteration ia to 
effectively greatly reduce the Reynolda number in the viacoua Jacobian 

matricea W and Zof equatlona (4.10b and c)[g]; thin haa no effect on the 

I 

Bteadyatate aolutlon, but can draatically affect the tranaient ao- 
lutlona^ The other ahortcoming of thia type of atart la that Intenae 
compreaalon and rarefacatlon wavea form, which, due to the Inexact 
nature of the boundary conditiona uaed in numerical aimulatlona, can 
^ cause long term distortion of the flow field [h]. | 

A second class of initial conditions is the same as the impulsive I 

start except that the velocities at solid bodies are gradually varied I 

from free-stream to the desired final values. This has been termed as 
a "penetration" initial condition, since, physically, this velocity 

I 

boundary condition implies that fluid must be sucked into and/or ejected | 

from the body. Initial conditions of this sort do not require a change j 
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In ctw chnraetnr of tho oot of oquatlono (4.10), howovor, thoy do r«- 
qulro that chay ba flaxlbla anough to hand! a auction and/or blowing 
froB tha boundarlaa. Siallar to tho Inpulalva atart eaaa, thara ara 
two aajor dloadvantagaa. Ona io that for a cloaad body ambaddad in tha 
flow fiald, ouch an an airfoil aaction, it ia axtranaly difficult to 
aatiafy global conaarvation lawa; thua tha aarly tranalant aolutiona, 
avan aftar tha daalrad boundary valocitiaa ara obtalnad, nay ba con- 
aidarably diffarant from axpactationa. Ralatad to tha firat potantial 
difficulty, tha othar problam with thin claaa of initial condition! is 
that if auction and/or blowing is not a natural boundary condition for 
a problam, than, again, tha transiant solutions for tho dasirod boun- 
dary conditions can not ba obtainod. 

Tha last typa of intial condition to bo considerod haro is what 
will ba callad a "gradual" start. This case is, in a sense, the re- 
verse of the Impulsive start in that the initial properties of the flow 
field and the boundaries are uniformly the same as the desired boundary 
condition on the solid body boundary. In most cases this means that the 
velocity field is everywhere zero. The start-up is achieved by 
uniformly accelerating the flow field everywhere except the "fixed" 
body boundary. I.i ,.£aotice, the easiest way to accomplish this is to 
add a body force term to the original equation set (3.1) and adjust 
the resulting difference equation set (4.10) accordingly. Details on 
how this modification is implemented are in Appendix I. 

C. Boundary Conditions 

Since the boundary conditions ultimately determine the form of the 
flow field, some care must be taken in formulating their numerical 
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•qulv«l«nt. This task la aspacially conplicatad by problama which hava 
no countarpart In tha analogous analytical fornulation. It is not cur- 
rently possibla to directly impose an "infinity" boundary condition on 
the numerical problem; thus, certain boundaries must be designated as 
Inflow and outflow boundaries and their boundary conditions determined 
so as to mimic the actual infinity boundary conditions. As mentioned 
previously (see discussion prior to equation (3.33)), an "extra" 
boundary condition is required on solid boundaries which must be deter- 
mined from the solution and not determine the solution. Also, re-entrant 
boundary conditions must be specified so that these computational boun- 
daries are effectively transparent to the flow field solution. 

Inflow boundary conditions are the easiest to formulate since they 
are normally placed far enough upstream of any solid bodies so that the 
flow in their vicinity may be assumed to be unaffected by their pres- 
ssnee. Thus, the inflow boundary conditions are completely determined 
by the corresponding infinity boundary conditions. 

Outflow boundary conditions, on the otherhand, are much harder to 
formulate, at least conceptually. Similar to the inflow boundary, this 
boundary is normally far downstream of any bodies, but, opposite to the 
inflow boundary, it must not affect the flow upstream of it. That is, 
this boundary must allow disturbances to pass "through" it without re- 
flecting them. One common method of at least partially accomplishing 
this is to determine these boundary values by extrapolation from the 
flow field solution, or, equivalently, by specifying a Neumann boundary 
condition. The technique chosen for this research was to use second 
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order extrapolation along the coordinate line which crosses the out- 
flow boundary. For example. If the Ith C-llne Is an outflow boundary, 
then 

- S-2.J 


q 


n+l 

I.j 




(4.11b) 


Re-entrant boundaries can arise either due to the analytical form- 
ulation (l.e., Couette flow) or due to the choice of coordinate system 
(l.e., polar coordinates). In either case a choice must be made be- 
tween treating these boundaries approximately, thus keeping the al- 
gorithm simple, or exactly. Typical of the approximate technique [b] 

Is to lag the values on this boundary for the solution of equation set 
(4.10) (henceforth tebroed Implicit boundary conditions), then to cal- 
culate these boundary values from an average of extrapolates of cor- 
responding points on the re-entrant boundary after the Interior field 
values are determined (henc forth termed explicit boundary conditions). 
In symbols this method could be Implemented as; 



Aq 


n-1 

I.j 


(implicit) (4.12a) 


Aq 


n 


I.J 




j+2 


2(Aq 


n 

I.j+1 


4 Aq 


n 

I,J+1 


) + Aq 


I .J+2-' 

(explicit) (4.12b) 


where ihe Ith C~line is a re-entrant boundary and j and J are the 
Indexes of the re-entrant n-line. 

The approach tak*n by this author was the more exact, but much more 
complicated, technique of continuing the algorithm across the re-entraut 
boundary. For cases in which the coordinate lines that crossed the 
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re-entrant boundary were periodic (l.e., polar coordinates), this only 
affected the inversion algorithm since the matrices were periodic 
block- triangular instead of Just block-tridiagonal. In all other 
cases, this "just*' means taking care that a block-trldiagonal matrix 
is formed in crossing the re-entrant boundary. In practice, the first, 
Inexact, method is preferable since it keeps the basic algorithm simple 
and much more problem independent. 

Body or wall boundary condi*:ion8 are the most difficult type to 
formulate, primarily because of their Importance, both analytically 
and numerically, in defining the flow. Unfortunately, this task is 
still mostly an art since a universally accepted method for treating 
the "missing" boundary condition on density does not exist. There is, 
of course, no problem in specifying the velocity boundary conditions 
in viscous flow; these are always assumed to be known quantities 
(usually identically zero). Also, the temperature boundary condition 
presents little trouble if it too is assumed to be specified (Isothermal 
being the easiest of these). 

It will be convenient to consider the implicit and explicit boun- 
dary conditions separately. For the implicit case, if the first 
n-llne is a body boundary, for example, then the trldlagnoal element 
-(Q^ 1 l^^^i 1 re-expressed in terms of known quantities 

and/or other differenced conservation vectors on the same ?-line. In 
general, this means determining the coefficient matrices H and K and 
the vector L in the linear combination: 

Aq" j = HAq" ^ + KAq" ^ + L (4.13) 
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which preserve* a block**trldlagonal system of equations when sub- 
stituted Into equation (A. 10c). These matrices were determined in 
this research by extrapolating the density along the C-llne and, in 
those cases In which the t«nperature was not explicitly specified, by 
lagging the temperature In time. The details of this technique are 
outlined In Appendix J; for other successful methods of evaluating 
these coefficient matrices, see Beam and Warming [c] or Steger [b]. 

After the solution of equation set (A. 10) the advanced boundary 
values on the body need to be determined. If a gradient temperature 
boundary condition Is specified by equation (3.32), then a straight- 
forward application of the appropriate difference approximations re- 
sults in the following periodic trl-dlagonal system of equation: 


-d) 


1.1 


1 - 1,1 


”l.l * ‘!>,>1,1 ■ "1.3 - '^"1,2 




(A.IA) 


To determine the surface density, the pressure gradient equation (3.33) 
was used. It also forms a slightly more complicated periodic trl-dlag- 
onal system: 


^Vi.! * <"r«"i,i * 


o ' j ** 1 + 1,1 “ ’^ 1,3 “ ^** 1,2 ^ ^2 


where the expansion of the tt^ and terms Is considered In Appendix C. 
For Illustrative purposes, the body has been assumed to be coincident 
with the first n-llne In both equations (A.IA) and (A. 15). Note that 
for bodies whose surfaces do not form continuous segments In the trans- 
formed coordinates, It was found to be advantageous to determine the 
value at the (re-entrant) end-points by the average of extrapolates 
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msthod outlined previously and then solving the resultant pure 
tri-diagonal system of equations for each segment. 


D. Artificial Viscoaity 

It is a well-known fact that, while the approximation of spatial 
derivatives by central differences, such as was done in developing e- 
quatlon (4.10), have many desirable properties, such as accuracy of 
derivative representation and simplicity of use, they make the al- 
gorithms containing them prone to nonlinear Instabilities. Thus, in 
any fluid problem in which a region of very rapid change occurs rel- 
ative to the corresponding coordinate density, some means of diffusing 
the resultant high frequency components of the solution must be ap- 
pended to the basic solution scheme. The technique used in this re- 
search was to alter the basic set of equations (4.10) by adding to the 
explicit equation (4.10a) the difference expression: 



+ (q 


1+2, J 
l.j+2 






(4.16a) 


everywhere except for points near boundaries where 

f ^‘*i+l,j ■ (4.16b) 

J 


is used. Note that the symbol q is used in the sense of equation (3.2). 
It can be readily shown that, provided the coordinate variation is 
"smooth", equation (4.16a) adds a term proportional to the fourth power 
of the local grid spacing and equation (4.16b) adds a term proportional 
to the square of this spacing. Thus, if e is proportional to At, for a 
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dens* enough mesh, these artificial diffusion terms should not seriously 
degrade the accuracy of the basic method except In those regions of high 
gradients (relative to grid spacing) where Its accuracy is suspect any- 
way. A linear stability analysis places an upper bound on the value of 
e of 1/12. Pulliam and Steger [g], from whom this dissipative technique 
was derived, also make use of an Implicit artificial viscosity appended 
to the W andQ matrices of equations (4.10b and c) ; however, this autlior 
found Its use to be required only for Impulsive type starts. 


E. Coordinate Derivatives 

There Is a serious potential problem with developing any numerical 
algorithm for solving the Navler-Stokes equations In conservative form. 
This can be descemed In equation (4.10a) by considering a region of 
flow In which the density, velocities, and energy can be taken to be 
constant. If proper consideration of the properties of the flux 
Jacobian matrices (equation (4.7)) is taken, then the first element of 
the vector RHs” Is proportional to; 


“ ^^n^l-l,J^‘^^^5^1,j+l ■ ^^5^1, j-1^^ 


Thus, since similar relations develop for the other elements of the 
vector RHS of equation (4.10a), if this algorithm is not to have spur- 
ious sources and sinks, the quantities in braces in expression (4.17) 
must be identically zero. In two dimensions, as is considered here, 
this can be assured by evaluating the coordinate derivatives (x^, x^, 
v . V ) by the same difference scheme used to difference the flux vec- 
tors A and B. As is pointed out in reference [b], this should be an 
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over-rldlng concern, even If the resulting differences are poor repre- 
sentations of the true derivatives (l.e., analytically generated 
coordinates should have their derivatives evaluated by the appropriate 
difference and not by direct differentiation). 
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V. COMPUTATIONAL RESULTS 


The types of flows modeled in this research can be divided into 
two clasaea on the basis of geometry. Couette flow, the first class 
of flows studied, was chosen primarily because it was the simpler of 
the two test cases that Beam and Warming [c] used to verify their 
original formulation of their approximate factorization algorithm. 

Thus, it was convenient to use this flow, both as a means of debug- 
ging the computer code form of the algorithm presented in the previous 
chapter, and as a means of verifying that the algorithm and computer 
code were consistent with the results of Beam and Warming. The other 
class of flow geometry, that of a NACA 64A010 airfoil section Imbedded 
in an other wise uniform stream of air, comprises the major thrust of 
this work. The particular airfoil section studied was selected on the 
basis of the number of experimental and numerical papers currently 
being published on transonic flow about NACA 64A010 airfoil sections 
at various angles of attack [b,J,k,l,m]. Also, the parameters of the 
flow were chosen to be representative of the cases studies in the ref- 
erences; thus, the Mach number (M-0.8), Reynolds number (Re*2,000,000) , 
Prandtl numbers (Pr«0.72 and Pr^-0.9), ratio of specific heats (>*1.4), 
and angles of attack (a*0‘’,2**,3.5‘*) were all set so as to facilitate 
comparison of results. In all of the cases discussed here, the three- 
point-backward version of the basic scheme was used (i.e. . 

A. Couette Flow 

Since the main purpose of this flow calculation was to reproduce 
the results of reference [c], the coordinate system and flow parameters 
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war* copied directly. The computational domain waa composed of a six 
by eleven (6 (-lines and 11 n-lines) rectangular grid with "movable" 
walla bounding the top and bottom surfaces (n"l and n*H) and with the 
left and right boundaries being re-entrant boundaries (i.e., (■! and 
are actually the same (-line). The flow parameters given by Beam end 
Warming are the Mach number (M-0.09), and the Reynolds number (Re>6.2) 
based on distance between the opposing walls; from their discussion. It 
was apparent that the working fluid used was atmospheric air, so the 
Prandtl number (Pr»0.72) and specific heat ratio (y* 1.A) were set at typ- 
ical values. 

In that the coordinate derivatives (( ,( ,n ,n ) are differenced by 

X y X y 

the algorithm used in this research Irregardless of the state of the 
fluid or boundaries, it was always felt to be Important that so-called 
"free-stream" trials be made in each case. For the Couette flow, this 
meant running a case with the fluid and walls initially in a uniformly 
quiescent state, and a case with both the fluid and the walls moving, 
in the plane of the walls, with an initial uniform velocity. In both 
cases the algorithm maintained the initial state for arbitrary lengths 
of time and for arbitrary At < 1. 

For the case of developing Couette flow, in which one wall is held 
fast, the fluid is initially at rest, and one wall is Impulsively brought 
up to the nondlmenslonal velocity of one . The basic test case was that 
for a time step of At=0.0116, which roughly corresponds to a Courant 
number of one. The coded approximate factorization scheme reproduced 
Beam and Warming's results exactly (possibly better than exactly) for 
this case, reproducing their conclusion that this method is excellent for 
determining transient solutions if the Courant number is less than unity. 
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Similar casaa vara run for Courant numbars of ten and one-hundred 
(At>0.08 and 0.8 reapectively) with excellent ateady-atate rcaulta and 
increaaingly poorer tranalent resulta - again the same experience ae re- 
ported in reference [c]. It ia thla author'a opinion that theae reaulta 
eatabllah at leaat the conalatency of thla formulation of the general 
approximate factorization acheme. 

B. NACA 64A010 Airfoil Section 

The coordinate aystem chosen for use with an airfoil with a rela- 
tively pointed trailing edge, such as this one, has proven to be of 
some importance in solving the numerical equivalent of the Navler- 
Stokes equations for the flow about the airfoil. There seems to be a 
consensus of opinion that the so-called "C" or "wake" coordinate sys- 
tem, shown schematically in Figure 1 , is most appropriate to a problem 
of this type [b,j]. Referlng to this figure, It Is seen that In this 
type of coordinate transformation, the lines of constant n form C's 
about the airfoil, with the inner-most line collapsing onto the airfoil 
surface and onto Itself to the right of the airfoil. While the lines of 
constant C extend from the outer n-line (n^J) to the inner n-llne 
(n*l), either terminating on the body or on the "cut" (the collapsed 
inner n~line) , If this coordinate system is to be useful, then care 
must be taken that the (;-lines which terminate on the cut pair-up with 
C-llnes on the other side of the cut so that they may be considered to 
be a single line running from the upper portion to the lower portion of 
the outer n-line. Looking now at the tranformed field, it is apparent 
that the left and right boundaries (labeled T* and F* will always 
be outflow boundaries, that the upper boundary (12) will be totally an 
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inflow boundary for O'* angla of attack and partially an outflow boun- 
dary otharwlaa, and that the portions of the lowar boundary labeled 
and Fj are re-entrant to each other. Although not indicated in the 
figure, thia coordinate ayatem formed a grid with 113 ^-llnea and 51 
n-linea; 20 of theae n-linea being concentrated within .05% to 2% of a 
chord length from either the airfoil or the cut. 

As diecuaaed in the Couette flow case, free-stream trials were felt 
to be an important test of both the solution algorithm and the coor- 
dinate system. For this coordinate system, runs were made fur flow 
field initial velocities of either zero or free-stream (including 
the body boundary, that is, the flow was forced to "penetrate" the air- 
foil). In both cases angles of attack were set at either O'* or 2'*, and 
the time step at At"0.01. Free-stream conditions were maintained for 
extended periods of time (t>l in all cases), within the truncation error 
limits, provided that the cautions of section IV. E were followed. That 
is, if the coordinate derivatives and the algorithm derivatives were not 
all of the same type (i.e., central differences), then free-stream, and 
Indeed a stable solution, could not be maintained. In particular, both 
tests of fourth-order differencing of the trailing edge coordinate de- 
rivatives and of the convective terms in the i ;<pllclt portion of the al- 
gorithm (see equation (4.10a)) proved to be destabilizing. 

The flow results obtained for the NACA 64A010 airfoil can be divided 
into several classifications based upon the method of starting the solu- 
tion process (impulsive, penetration, or gradual), the turbulence model 
(laminar or algebraic Baldwin and Lomax model) , and the angle of attack 
(0°, 2°, or 3.5"*). Initial cases were all run as laminar flows at zero 
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d«gr««s angle of attack, thua only being dlatingulahed by the type of 
Btart utilized. Later caeea were all etarted and/or changed gradually 
and employed the algebraic turbulence model after an initial laminar 
Btart ing period; thua theae caaea are dietinguiehed only by the angle 
of attack chosen. 

One of the first starting techniques tried was the impulsive start 
method. Similar to the results of Steger and Bailey [n], numerous 
trials established that addition of only explicit artificial viscosity 
to the numerical algorithm was not sufficient to overcome the non- 
linear instabilities produced by the sharp gradients generated by this 
type of start. Although Steger and Bailey were able to successfully 
utilize this technique by also including implicit artificial viscosity, 
it was the experience of this author that no combination of explicit 
and implicit artificial viscosities would lead to suppression of the In- 
stabilities inherent In this method of starting the fluid flow. It is 
conjectured that this is due to the differences in the meshes in either 
case; although both were "C" type coordinates, the ones used In this 
work had much smaller cell sizes near the leading edge. Thus the 
courant number was appreciably larger in the leading edge region In this 
research as compared to that of Steger and Bailey. If this conclusion 
Is correct, then use of a smaller Initial time-step, At , should allev- 
iate the problem. However, this was not attempted since a primary goal 
of this project was to develop techniques which allow the use of larger 
time-steps. 

Figure 2 Illustrates the transient behavior of the flow field for 
this type of start. Notice that even for this late a time (0.6 non- 
dimensional time units) that the flow, outside of a very thin region 
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n«4ir the airfoil and tha waka cantar-lina, ia eaaantially undiaturbad. 

Tha axcaptiona being the wake expanaion at 0.6 of a chord from the 
trailing edge and tha coa^raxaion wave anananating from tha leading edge. 
Alao note the onaet of fatal nonlinear inatability indicated ty tha 
"wigglea" in the Mach 0.75 contour denoting the compraaaion wave. 

A more aucceaaful starting technique was tha penetration start. 

As outlined earlier in thia paper, this start was accomplished by using 
a fifth-order poly.iou'ial to reduce the suction and blowing at the air- 
foil surface from free-stream values to no-slip values in one nondimen- 
sional time unit (one-hundred time-steps at the standard time-step 
used throughout this work of At"0.01). This method was fully success- 
ful in getting the flow started, however, as noted previously, it had 
the disadvantages of requiring non-slmple body boundary conditions and 
of causing early boundary layer separation over the aft portion of the 
airfoil. Ultimately, this approach was abandoned since it was felt that 
the last mentioned deficiency of this method could lead to a valid, but 
undeslred, steady-state solution. That is, since the fluid flow Itself 
can take on a variety of steady flow states for Identical boundary con- 
ditions depending solely on the Initial conditions [o], it was felt 
chat the early, forced separation of the boundary layer due to this 
technique could lead to a steady-state not found under normal flight con- 
ditions. Also since there were more grid points on the forward section 
of the airfoil than there were on the aft section, the lack of global 
conservation of flow properties alluded to in a previous chapter pro- 
duced by this starting technique lead to extreme cooling of the trailing 
edge region which took an undesirable number of time-steps to relax. 
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Pitur* 3 and Figure 4 present the Mech number contours of the 
trenelent flow developed midway through end at the end of the penetration 
■tart, respectively. Coshering these figures with the Impulsive start 
Mach contours In Figure 2. It Is Inmedlately apparent that the pene- 
tration start allows a much larger region of the flow field to devel- 
op for the same time Increment. Also note that the various expected 
flow regions are already well-developed by the end of the start-up 
period. However. It Is also well to notice that the forward deceler- 
ated flow region and the surface accelerated regions are both distorted 
In the downstream direction; and that the wake contains an appreciable 
highly decelerated region In its core. 

In that the single laminar flow case treated by this author was 
initiated by the penetration start method, the gradual start technique 
will be discussed following the presentation of the laminar flow re- 
sults. These results are summarized in the Mach contour plots, Figures 
4 through 9. Several Interesting observations can be made by an In- 
spection of this flow's time history. First note how the retarded flow 
region expands from its downstream skew due to the penetration start, 
to an overshoot upstream skew, and then contracts to the expected, 
slightly elliptical contours. Looking at the aft boundary layer and 
near wake, it should be noted that these regions first thin out from the 
thickness induced by the blowing during the start-up; then, between non- 
dimensional times two and three, separation occurs resulting In a pro- 
gressively thicker separated region and near wake region. Coincident 
with the occurance of separation, the accelerated regions start to 
spllt-up Into two zones of maximum Mach numbers. Within each zone, the 
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Mach nuBlara dacraasa with tiua, although tha portion of the flow within 
each sone remalna relatively conatant aa the outer Mach contour a 
"pinch-off". The forward sone prograaaae toward the leading edge aa the 
aft zone llfta-off with tha eeparated boundary layer. 

The cccurance of theaa later phenomena can be explained ae logical 
conaequencea of the atrongly eeparated flow. The flow field outalde of 
the boundary layer and eeparated region "aeea" the airfoil aa if the 
eeparated region were a continuation of the airfoil aurface; thua the 
outer flow muat flrat accelerate around the leading edge region and then 
over the eeparated region. Similarly, the general drop in maximum Mach 
number, from a high of 1.15 before separation occured to 1.05 in Figure 
9, and the contraction of the Mach contoura delineating the accelerated 
region, can be attributed to the eeparated boundary layer giving the 
airfoil the appearance of a thin wedge to the outei flow. Unfortu- 
nately, this author was unable to obtain any comparative studies for 
this case (not to say that some might not exist). 

The gradual starting method in all ways turned out to be an almost 
perfect flow initialization technique. It only required a modest amount 
of explicit artificial viscosity, smoothly accelerated the airfoil 
at 0.0 and 3.5 degrees angle of attack from rest to Mach 0.8 within a 
nondimensional time of one, and even provided the bonus capability of 
changing the angle of attack in mid-solution. Best of all, in this 
author’s opinion, it al.jws the most realistic starting procedure of all 
those considered in this research (including some not actually utilized, 
such as potential flow solutions). In all applications the same fifth- 
order polynomial as that used in the penetration start to vary the 
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■Irfoll surfaca velocltlaa waa uaed to vary the free-atream velocities 
from zero ^o the deaired ateady-atate valuea (this was handy and also 
allows some comparison of the penetration and gradual starts to be made) . 

The first application of this method was to the NACA 64A010 air- 
foil at zero angle of attack, as illustrated by Figures 10 and 11. A 
comparison of the three starting techniques, Impulsive, penetration, and 
gradual (Figures 2, 3, and 10) Is Instructive, but note that the free- 
stream Mach number of the flow plotted In Figure 10 Is 0.73 rather than 
the Mach number of 0.8 for the other two. Contrasting the Mach con- 
tours of Figures 4 and 11, which Indicate the floi^ fields at the end of 
the penetration and gradual starts, respectively, one sees that the 
flows are very similar except In the aft boundary layer and near wake 
regions. Here the gradual start produces a much thinner boundary layer 
and near wake, except for the Mach 0.75 contour which has the bubble 
shape characteristic of experimental results [m]. 

Immediately after the completion of the gradual start (i.e., at 
nondlmenslonal time of one) , the Baldwin and Lomax turbulence model was 
turned on and the solution continued out to a non-dimensional time of 
six, as with the laminar case. A comparison of the turbulent results 
(Figures 12 through 16) with the laminar results (Figures 5 through 9) 
proves to be enlightening. First of all, note that the forward decel- 
erated region develops almost Identically in both cases, which is as 
would be expected. Also, the location and Intlal slopes of the upstream 
accelerated flow Mach contours 0.85, 0.90, 0.95, and 1.00 remain in good 
agreement between the two flows until after the nondimens ional time of 
four. Even after this point, the first three contours agree in upstream 
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location and initial slope near the airfoil as far as both cases were 
Cakcnt This particular behavior of both the decelerated and accelerated 
regions can be attributed to the fact that the turbulence model, as lm> 
plemented by this author, does not become effective until slightly down- 
stream of the first minimum surface pressure location (l.e., turbulent 
transition is modeled), thus, upstream of this point one would expect 
laminar and turbulent flows to behave similarly at least for early 
times and at the freestream Mach number chosen (0.8). A final obvious 
contrast is the presence and non-appearance of separated flow in the two 
cases, which reflects the elementary experimental and theoretical pre- 
cept that turbulent flows can tolerate an adverse pressure gradient 
much better than laminar flows. 

In comparing the turbulent Mach contour results It is first of 11 
apparent that the near wake reaches a quasi-steady state very early, the 
only subsequent significant changes being the thickening and ultimate 
closure of the middle wake ret,lon. This last result is somewhat puz- 
zling, especially since the entire wake shows practically no change at 
all between nondlmensional times of four and five. Also note how the 
accelerated flow Mach contours progressively expand out from the airfoil 
and gradually shift toward the leading edge. Their initial downstream 
skew also tends to turn into an almost symmetrical appearance by the 
last time step plotted (Figure 16) . The Mach one contour never be- 
comes normal to either the chord line or the airfoil surface and the 
spacing of the downstream Mach contours never take on a crowded appear- 
ance. Thus it is hardly justified to claim that a shock develops on 
the airfoil; however, it appears from the Interferogram of reference [m] 
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that in the experimental case the shock la rather weak. None-the-less, 
It was decided at this point that the ability of the numerical algorithm 
to correctly develop shocks needed to be given a less ambiguous test. 

Referring to either reference [k] or [m], It Is clear that in ex- 
perimental studies of this airfoil, a moderately strong shock appears 
on the upper surface for an angle of attack of two degrees. Thus, this 
angle of attack was chosen for the next test case which used the pre- 
vious results at nondlmenslonal time of four for Initial conditions. 

The gradual start method was used to change the angle of attack from 
zero to two degrees, effectively by causing the airfoil to "fall" at a 
speed sufficient to produce this change (note that this is not the same 
as rotating the airfoil, and thus, will produce different transient 
states and could produce a different steady state) . Then the solution 
was continued out to a nondlmenslonal time of seven, at which point the 
solution process was terminated due to the financial suicide of the 
university's computer center. 

Figures 17, 18, and 19 contain Mach contour plots similar in nature 
to those discussed previously. It is immediately noticeable that the 
flow in all regions is no longer symmetrical about the chord line. The 
progressive changes in the decelerated flow region and the wake region 
are as in the zero degree angle of attack case except for their respec- 
tive skew in the upstream and downstream directions. More notably dif- 
ferent is the evidence of much more strongly accelerated flow along the 
upper surface; however the form and density of the Mach contours still 
don't suggest that a shock is present. To clarify this lack of an ex- 
pected flow phenomina, the Mach contours plot Identified as Figure 20 
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was made with the range of Mach numbera being 0.96 to 1.04 and the In- 
crement being 0.01 rather than the range and Increment uaed In the other 
Mach contour plots. Aa can be seen, there Is a pronounced difference 
between the upper and lower contour groups, as expected. However, the 
angle and density of the upper groups downstream foot are not as would 
be expected for even a rather weak shock. Although; as mentioned pre- 
viously, It was not possible to continue this solution so as to check 
this hypothesis. It was conjectured that much of the spread In these 
Mach contours was due to the relatively large amount of explicit art- 
ificial viscosity used. This was reasoned from the fact that the so- 
lution showed no signs of "wiggles", which is unusual for numerical 
solutions of flows containing shocks. 

A valuable check of the internal consistency of the numerical 
solution, as well as an additional means of comparing the results with 
the available experimental data, can be obtained from a density contour 
plot such as Figure 21. First of all, comparison of these contours 
with those of reference [k] reveals at least good qualitative agreement 
with the "set" angle of attack of two degrees data set. Although ref- 
erence [m] does not present results for this particular angle of attack, 
the density contours of Figure 21 do seem to represent a case between 
the angles of attack of G and 3.5° that are shown. These comments must 
be tempered with the fact that both references clearly show a shock on 
the upper surface of the airfoil whereas the Mach number and density con- 
tour plots of the algorithm used by this researcher clearly do not show 
such a shock. 

Some Interesting observations can be made In comparing the density 
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contours of Figure 21 with the Mach contours of Figure 19 which liave 
some bearing on the Interpretation of experimentally determined density 
contours obtained by Interferometer, such as those of reference [m]. If 
one makes the same assumption as used In the cited reference, that the 
flow Is Isentroplc everywhere except In the boundary layer and across 
any shocks, then one can obtain the following relationship between non- 
dimensional density and Mach number: 

- 1 + ^(m2 - m2) (5.1) 

where for this equation only, M is the local Mach number and M^ is the 
reference Mach number (0.8 In this case). (For the convenience of the 
reader, a rough rule of thumb summarizing this equation for the current 
application is that a change In density of 0.05 - one density contour 
of Figure 21 - results In a change of Mach number of about 0.06 to 
0.08). Using this relation. It Is easy to see that there Is good agree- 
ment between the two sets of contours except for four flow regions. The 
most puzzling discrepancy Is the lack of alignment of the Mach 0.75 and 
density 1.05 contours upstream of the leading edge, since this region 
would seem to be nearly Isentroplc. Possibly this Is due to the fact 
that the flow has not quite reached steady state. Similarly, the lack 
of agreement between the Innermost accelerated region contours of the 
two plots could be attributed to a combination of nearness to the non- 
Isentroplc boundary layer and the transient nature of the flow; however, 
these arguments are not convincing. More importantly, note how the 
downstream intersection of the accelerated flow region density contours 
with the boundary layer are shifted upstream and are more nearly normal 
to the surface than the corresponding Mach number contours. Also note 
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the almost total lack of similarity of the Mach 0.75 and density 1.05 
contours near (or In) the wake - clearly due to the non-laentroplc 
nature of the wake region. 

The last case considered was that of the same NACA 6AA010 airfoil 
at an angle of attack of three and one-half degrees. This case was con- 
sidered due to Its moderately strong shock development and because of 
the excellent set of wind tunnel data gathered for this configuration 
by Johnson and Bachalo [m]. The results for this case are In many ways 
the most successful of those attempted In that a distinct shock Is seen 
to form; however, there exist a number of unresolved problems connected 
with this solution. Most disappointing of which was the unexpected on- 
set of numerical Instability around a nondimension time of six. 

Figures 22 through 27 Illustrate the Mach number contours generated 
at various time Increments by this solution. Again, the gradual start 
method was used, completing at nondlmenslonal time of one. In this 
case, though, the airfoil Is Initially at an angle of attack. Comparing 
the progression of the Mach contours, one is struck by a similarity be- 
tween all of the cases studies in this work; that is that all of the 
flow regions tend to expand from the start until sometime between non- 
dimensional times four and five, at which point they start to contract, 
except for possibly the accelerated region on the shock side of the air- 
foil. The reason for this surely has something to do with the rapidity 
of the start; it is doubtful that the outer boundaries contribute to 
this phenomlna since they are ten chord lengths away from the airfoil 
leading edge in all directions. 

Since much of the description of the results for this case para- 
llels that of the previously presented cases, this presentation will 
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concentrate on the two distinctive results of this solution. The most 


heartening of these Is the development of what can only be described as 
a shock on the upper surface of the airfoil. The Mach contours on the 
upper rear portion of the airfoil section are clearly crowding together 
Into a thin bundle which Is normal at least to the chord line. This 
shock formation can also be observed from the coefficient of pressure 
plots (Figures 28 through 33) In that the downstream portion of the up- 
per surface curve progressively steepens. Lastly, the obvious signs of 
boundary layer separation on the upper airfoil surface near mid-chord, 
also gives support to the development of a shock-llke pressure gradient. 
The other distinctive feature peculiar to this solution Is the appear- 
ance of the wake "bubble" contour just off the trailing edge. While no 
definitive explanation of this effect can be put forward, It Is specu- 
lated that this anomaly Is either due to the onset of boundary layer 
separation, a "bug" In the coded algorithm, or the unexpected, but ac- 
tual way a flow of this type would physically occur. Unfortunately, as 
a comparison of the Mach contours and coefficient of pressure plots 
with each other and with the experimental results of reference [m] in- 
dicate (see Figure 33), this flow Is approaching, but has not yet 
reached, steady state. However, due to the aforementioned stability 
problem and the deplection of the time allotted to this research, this 
work must end on a promising, yet frustrating note. 

As an afterword to those who may wish, tor whatever teason, to re- 
produce or continue this work: in addition to the fixed flow param- 
eters mentioned at the first of this chapter, all cases were attempted 
under the adiabatic wall assumption, with the non-dimensional time-step 
always fixed at At«0.01, with the explicit artificial viscosity 
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coefficient set at e*8At except for the last two-hundred time-steps of 
the last study where c"4At was used, and with the coordinate system 
generated by Dr. Joe Thompson, Aerospace Department, Mississippi State 
University, Missiasippi 39762. 


48 


VI. CONCLUSIONS AND RECOMMENDATIONS 


On the whole, the primary objective of this research work was 
obtained, although some subsidiary goals must be left to future inves- 
tigations. An approximate factorization algorithm modeling the full, 
Reynolds-averaged, Navier-Stokes equations was successfully developed 
and Implemented in a general-purpose computer code. This computer pro- 
gram Includes the desired features of an algebraic turbulence model, 
penetration and gradual start-up routines, variable Inflow option, 
variable body surface velocity capability, body surface heat flux or 
temperature distribution specifications, implicit and/or explicit 
artificial viscosity options, and translational rigid body motion 
capability. In addition, the basic algorithm was formulated in con- 
servation form on a body-fitted curvilinear coordinate system. 

The amount of computer time required per time step averaged 
35 seconds on the Univac 1100/80, and 4.5 seconds on the CDC Cyber 203 
(this works out to approximately 8 hours and 1 hour total computer 
time for a complete time-history of the flows considered here) . How- 
ever, although the cases studied are indicative of the full capability 
of the algorithm as conceived by this researcher, further work is 
clearly required in the area of model verification. 

Of the multitude of observations that can be drawn from the body 
of this research, it is felt that the following are most worthy of 
attention. While the body-fitted coordinate technique was fundament- 
al to the Implementation of this algorithm, two areas in which im- 
provements are urgently needed can be identified. These involve the 
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development of rational and efficient means for determining coord- 
inate line spacing and placement and for dynamically concentrating 
coordinate lines in the region of developing shock waves. As nested 
earlier in the text, and related to the last topic, the standing prob- 
lem of how to express truncation errors in curvilinear coordinates in 
as concise and meaningful a manner as the "order of" method of Car- 
tesian coordinates, should be resolved if differencing in transformed 
coordinates is ever to have a rational basis. Boundary condition 
problems tended to dominate the difficulties encountered in application 
of this algorithm. In particular, the much-neglected downstream 
boundary condition requires more thought than has been given it to 
date - in fact, the entire problem of replacing "infinity" boundary 
conditions by finite distance boundary conditions needs a firmer theo- 
retical basis. In using the adiabatic boundary condition in this re- 
search, numerical evidence seems to indicate that, at least in the im- 
plicit portion of the algorithm, either the density or the temperature 
can be extrapolated into the field, but not both together if nonlinear 
InstabilltleQ are to be avoided. This points to the whole unsatisfac- 
tory situation in which there is no firmly established "fourth" body 
boundary condition - a topic which needs some more fundamental re- 
search efforts. Lastly, the results of this research indicate that 
construction of flach number contours from experimentally determined 
density contours may not always accurately reflect the true Mach number 
contours. Thus, it is suggested that published experimental results 
Include the density contour tracings in place of, or along with, the 
calculated Mach number contour plots. 
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b. Transformed Field 


Figure 1. Wake Type Coordinate System 
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MACH NUMBER 



Figure 2. Mach Number Contours - Impulsive Start 
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Figure 4. Mach Number Contours - Penetration Start 
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Figure 5. Mach Number Contours - Penetration Start 



Figure 6. Mach Number Contours - Penetration Start 
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Figure 6. Mach Number Contours - Penetration Start 


MACH NUMBER CHi^iGES 
BY 0.05 BETWEEN CONTOURS 
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Figure 9. Mach Number Contours - Penetration Start 
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Mach Number Contours - Gradual Start 
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Figure 11. Mach Number Contours - Gradual Start 
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Figure 12. Mach Number Contours - Gradual Start 
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Figure 13. Mach Number Contours - Gradual Start 
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Figure 14. Mach Number Contours - Gradual Start 
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Figure 16. Mach Number Contours - Gradual Start 
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Figure 17. Mach Number Contours - Gradual Start 



MACH NUMBER CHANGES 
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Mach Number Contours — Gradual Start 
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Figure 19. Mach Number Contours - Gradual Start 



Figure 20. Mach Number Contours - Gradual Start 
ot = 2°, t = 7.0 - Expanded Scale 
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Figure 21. Density Contours - Gradual Start 


MACH NUMBER CHANGES 
BY 0.05 BETWEEN CONTOURS 



Figure 22. Mach Number Contours - Gradual Start 
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Figure 23. Mach Number Contours - Gradual Start 


MACH NUMBER CHANGES 
BY 0.05 BETWEEN CONTOURS 



Figure 24. Mach Number Contours - Gradtial Start 



Figure 25. Mach Number Contours - Gradual Start 




Figure 27. Mach Nus&er Contours - Gradual Start 










APPENDIX A 


Nondlmmiionallzat ion 


Although the uae of nondlmenalonal equationa haa many wall-known 
advantages. It has the major draw-back that there Is no unique dimen- 
sionless form for any given set of equations. Thus, two otherwise 
formally identical sets of equations may appear to be different solely 
due to the nondimenslonallzation used. Slrce, in this author's exper- 
ience, the only consistent criteria for nondimenslonallzlng the 
Navier-Stokes equations is that the Reynolds number (Re) and Prandtl 
number (Pr) should appear in certain traditional locations, this has 
been taken as the primary objective of the nondimenslonallzation used 
in this work. The advantage of this criteria is that it keeps the set 
of partial differential equations from becoming cluttered; on the other 
hand, the suppressed dimensionless ratios now show up in the auxiliary 
algebraic relations and the boundary conditions. 

All of the dimensionless variables and ratios ultimately are ref- 
erenced relative to the following characteristic parameters: 


Parameter 


Name 

Length 

Velocity 

Density 

Shear Viscosity 
Thermal Conductivity 


Dimensions 

L 

L/T 

M/L^ 

M/L-T 

M-L/e-T^ 


Specific Heat at Constant L^/9-T 
Pressure 

Specific Heat at Constant L^/6-T 
Volume 
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Paranef r 


Name 


Dlmenalona 


Speed of Sound L/T 

where the dlmenalona are length (L), time (T), maaa (M), and temp- 
erature (6). The apeed of aound, C , la uaed only In altuatlona 

o 

where the Mach number would be more natucal than a dlmenalonleas temp- 
erature. Theae reference ^uantltlea are then used to form the following 
Hat of nondlmenalonal variables or constants: 


Kondlmenslonal Variable 

u,v 

p 

M,X 

k 

t 

P.e 

T 

q" 


Reference Quantity 
I 


1/u 


p u^ 

O 0 

u2/C 
o p 

3 

p U^ 

O O 


with the Sutherland reference temperature, , appended to the charac- 


teristic parameter list, the following nondimensional characteristic 
flow pi^.rameters can be formed: 

Reynolds number 
Prandtl number 
Mach number 

Ratio of specific heats 
Sutherland reference temperature 


Re - u /pu 1 
o o 

Pr » C p /k 
p o o 


M 

Y 


u /C 
o o 

C /C 
P V 


S, « CpS /u^ 
1 o o 
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APPENDIX B 


Coordlnaf Transformation 


It Is widely accepted that any finite difference approximation to 
the Navier-Stokea equations should retain the global conservation prop- 
erties of the Integral form of these equations. This feature can be ob- 
tained, for usual differencing, only If the transformed equations retain 
the same general form as the differential Navler-Stokes equations in 
Cartesian coordinates. In this paper, this was accomplished In trans- 
forming equation (3.1) Into equation (3.15) by use of the following re- 
lation (note the Einstein summation convention Is used) : 


9A£ dAl 

where _1 3£i . 

1 “ J 

and J Is the Jacobian of the coordinate transformation, 
plication under consideration, the vector components, A^ 
ordinates, and , can be identified as: 


(B.l) 

(B.2) 

For the ap- 
, and the co- 


X, - t 


^1 “ ^ 


A 2 - A + D 
Aj “ B + E 


X 2 ■ X 


5 

n 


(B.3) 


and the Jacobian, J , is: 


J = 


1 - T T 
t X y 

^t ^x S 


n n 
X y 


X 


C n - C n 

X y y X 


(B.4) 


since Tit. 
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Now relation (B.l) can ba aatabllahed by first drveloping the 
following chain of equalities: 


-c 

3x 


C, - C 


j ''4k '’4 4k 9x^ '’i '^4 


J6. 


i4 3C^3x^ ^4 " 

where is the cofactor of the 4,k element of the Jacobian determinant 
(B.A) and the first equality ia proved in most standard tensor treatments 
(e.g. Sokolnlkoff [l], page 103). Then expansion of the right-hand-side 
of equation (B.l) gives: 

^ 35^ ^i ^ 3C^^J ^j3x^ ^i^ ’ J 3Xj 3^3Xj ^1 3Xj 

(B.6) 

which on use of the last equality of (B.5) reproduces equation (B.l). 


(B.5) 
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APPENDIX C 


Normal Pre««ure Derivative 

The particular form of Che normal preaaure derivative equation 
used In this work (equation (3.33)) can most easily be developed from 
the Navlcr-Stokes equations In curvilinear coordinates, equations 
(3.21) through (3.26). First the x- and y- momentum equations are 
expanded to: 

u[pJ+(p*U)^ + (p*V)^] + p*[u^+Uu^+Vu^] + - Rj (C.l) 

v[pJ+(p*U)^ + (p*V)^] + P*[v^+Uv^+Vv^] + j[5yP^+nyP^] - R 2 (C.2) 

where p*** •• p/J and Rp R 2 are the unexpanded viscous terms of the x- 

and y- momentum equations respectively. Now the first term In both 

equations (C.l) and (C.2)vanlshes due Co the continuity equation; thus. 

If equation (C.l) Is multiplied through by n and equation (C.2) by n 

X y 

and then added, the result is: 

■"x'l ^"2 

Finally, If this equation is multiplied by 3! /a and the equation of 
state, P ■ (y“1)pT/y , is used, then equation (3.33) results. 

The finite difference form of equation (3.33), the normal pres- 
sure derivative equation, used in this research, equation (4.15), was 
put into the particular form chosen so as to facilitate the use of a 
common solution algorithm for pressure and temperature when the normal 
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prasBura derlvatlva la Bpeclflad (aaa aquation (4.14)). Tha ir terma, 
avaluatad at 1,1 ara: 

”2 ’ ’ 

Tha derlvatlvea are approximated by the appropriate finite difference ex- 
preaalon from Appendix G and evaluated prior to the solution of equation 
(4.15). 
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APPENDIX D 


Boundary Fltt«d Coordlnatos 


Since all of the problcma conaldered In thla paper poaaeas single 
body, two-dimensional geometrios, and since this class of coordinate 
generation problem has been extensively examined in the literature (see 
[d,e,h,j], for example), this appendix will be devoted to clarification 
of equations (3.37) and (3.38) and a summary of the coordinate contrac- 
tion technique discussed more completely in reference [j]. The symbols 
used in the aforementioned equations are defined to be: 


K * y 
n 'n 


(D.l) 


®c • Vn* 


(D.2) 


Y ■ X^ + 


(D.3) 


J - J ■ X y - X y 
c 5 n n C 


(D.4) 


The boundary conditions on the coordinate generation equations 
(3.37) and (3.38) are determined by specifying the ^ and n distribution 
on the body and the inflow and outflow boundaries. As mentioned in the 
text, normally this is done so that each boundary is a line of constant 
^ or n. The system of partial differential equations is then solved by 
approximating all derivatives by "second order" central differences (see 
Appendix G) and using a point successive over-relaxation iteration scheme 
to determine the solution of the resulting nonlinear simultaneous dif- 
ference equations. 

Of course, before the coordinate generation equations can be solved, 
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th« coordinate lino attraction functions, P and Q , mist ba specif led. 

c c 

The general effect that these functions can have on the reaultlng co- 
ordinate systm, and cautions about the forms they can take are fully 
developed in reference Cd]; it is enough to note here that these at- 
traction functiona can have a dramatic effect on the coordinate spacing 
near a boundary. Since the high Reynolds number flows considered in 
this work are very dependent on boundary layer interactions, it was de- 
sired that the coordinate spacing be very denee inside the boundary 
layers and, thus, near the body boundary. As noted in references [j,k], 
this can be readily accomplished by taking the appropriate attraction 
functions (Q^ if the body surface is a line of constant n) as: 


Qg(C,n) 


(2-fninK) 

(l+ninK) 


inK 


0.5) 


where the spacing parariit T K is determined from: 

NK^“^ - 1 in^tT^ 


- 1 




(D.6) 


The parameters in this last expression are defined so that If the body 
surface coincides with the first n-line, then r^ is the radius of a 
circle which circumscribes it, and if the outer boundary lies on the Jth 
n-line, then r^ is the radius of a circle which Is tangent to It. This 
set of equations has the effect that the Nth n^line will be approximate- 
ly a distance from the body surface when they are used in the co- 

ordinate generation scheme. In this study, the 2nd n-line was specified 
to be one-hundredth of the boundary layer thickness from the body. So, 
by approximating the boundary layer thickness from the laminar Blaslus 


flat -plat* boundary layar aolution, thl* la aqulvalant to 


fjj • O.Ol 




(D.7) 
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APPENDIX E 


Plux V«ctor» and Th«lr Jacobian Mbt.rlc« » 


Since ch« dynamic flux vactora (A and B) and Cha viacoua flux vac- 
tora (DpD 2 ,R| and E 2 ) ara ainilar in form within aach class » 1st the 
vectors A and D ba represantativa of each class in the following dis- 
cussion. If for the present, the transformed coordinates (C and n) are 
taken to be constant in time, than the chain-rule nay be used to form 
the following identities: 

- Pq^ (E.l) 

- itq^ + Wq^^ (E.2) 


where the Jacobian matrices P, n, and W are defined to be 



TT 


()D ,, 3D 

• W - -r— - 
?q 3q 


Now the second identity (E.2) can be put in the form: 


(E.3) 


- (n-W^)q^ + (Wq^)^ (E.4) 

If the viscous coefficient (i.e., u times expressions (3.10) and 
(3.11)) are taken to be locally constant in time, which ultimately intro 
cudes an error proportional to their first time derivative times At^, 
then it can be shown by direct evaluation that and from equation 

(E.4): 

Equations (4.5) then directly follow from 

AA - Aj.At " Pq^At = PAq + O(At^) (E.6) 

AD = D^At = (Wq^)^ = (WAq)^ + 0(At2,p^At) (E.7) 
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Under the e..y>unption used In ••tabXishing equation (E.S), 
equation aeta (4.6) and (4.7) can alao ba confirmed by direct evalu- 
ation. It ia intareating to nota that for any flow in which tha via- 
coaity ia aolaly a function of twaparaturc (auch aa laminar flow with 
a Sutherland viacoaity law) • than tha flux vactora are homogenaoua 
functiona of varioua dagraaa in tha conaarvation vector, q , and ita 
firat derivativea; thua equationa (4.6) and other intaraating ralation- 
ahipa can ba aanily ahown. 

Aa a final nota, if the tranaformed coordinatea are known functiona 

of time, or at laaat are advanced in tima prior to aolving for the flow 

variablea, then it can be ahown that if the derivativea 

n ,n ) are evaluated at time (n4'e)At, where 0 ia aa introduced in e- 
X y 

quation (4.1), then thaaa derivativea can be treated aa if they were 
locally conatant with time without degrading the accuracy of the ap- 
proximate factorization scheme. 
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APPENDIX P 


Explicit Form of th« Jacobian Matrlf 


These matrices are found by a straight-forward, If rather tedious, 
process of re-cxpresslng the vectors of equations (3.22) through (3.26) 
explicitly In terms of the conservation variables (q^ p/J, q 2 ■ pu/J, 
q^ * pv/J, and q^ ■ e/J, and then preforming the differentiations In- 
dicated by equations (4.5). They are: 


Q - 


W 


0 


X 

^y 

0 


c ♦- 

X 

uU 

U-(y-2)C^u 

tyU-(Y-l)CjjV 

(Y-I)^., 


K ♦- 
y 

vU 

5jjV-(Y-l)CyU 

U-(Y-2)tyV 

(Y-l)ty 


(♦-T 


5^Ts-(y-DuU 

CyV(y-i)vU 

yu 

4SI 


■ 

0 



'’y 

Si 

0 


n 

X 

uV 

V-(Y-2)njjU 

n u-(Y-l)n V 

y * 

(Y-l)n, 


n 't'- 

y 

vV 

njjV-(Y-l)nyU 

V-(Y-2)riyV 

(Y-DOy 


(♦-T 

m 


n T^-(y-1)uV 

X 8 

nyV(Y-l)vV 

yV 

mm 



‘o 


0 0 


o" 

1 

-(a 

jU+a2v) 

®2 


0 

P 

-(a2U+a^v) 

^2 ®3 


0 


-(W 


-(Wjj+a^v) 





0 0 


o’ 

J. 

-(bjU+b2v) 

^ b2 


0 

P 

-(b 

ju+b^v) 

b3 


0 


-(X 

^2U+X43V+b5E) 

-(Y2i+b5u) -(Y3^+b5u) 



(F.l) 


(F.2) 


(F.3) 


(F.4) 


Y - - 
P 


Z - 


-(bjU+bjV) 


'41 


0 

b, 


-(b2U+b^v) b^ 


-(CjU+C2v) 

-(C2u+Cjv) 


0 

b. 


-(X2j+b5u) -(Xjj+bjV) bj 




0 

0 

0 

c. 


(F.5) 


(F.6) 


where the contravarlant velocities, U and V, are given by equations 
(3.28) and (3.29) and the viscous coefficients are the same as those de- 
fined by equations (3.3C). The subscripted Jacobian matrix symbols in- 
dicate the appropriate element In their explicit expansion (e.g., W 2 ^ ■ 
-(a,u+a„v)). The other symbols introduced here (4>, T and E) are de- 
fined to be: 

4> ■ • i) (u^+y2) 

T - yE-4' (F.7) 

8 

E “ e/p 
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APPENDIX G 


Difference Approximations 

The finite difference expreeslone used In this work are pre- 
sented below. In all cases the functions a(^,n) and f(^,n) are duorniy 
functions representing typical actual derivative expressions encount- 
ered In the development of the text. Equations which have an obvious 
analog when the variable of differentiation or the order of differen- 
tiation la changed have not been repeated. 


Central Differences: 


1 


2^^1+r^l-l^ 


*55 ■ 


'5n ■ 

4^^i+i ,j+r^i+i ,j-r^i-i ,j+i"^^i-i , j- 

<“'5>5 • 


'“'n>5 ■ 

4^"i+i,j^^i+i,j+r^i+i,j-P'“i-i,j^^ 


(G.l) 


One-Sided Differences: 




‘55 

n n 


-<^-3-^'l-2^5f^.j-2f^) 






-Uj (f j^^-4f j^2+5f j^l-2f j ) 
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n 5 




“ 2^ ^“i+1 , J" Vl , J ^ ^^1 , j+2"^^i ,J+l‘''^^i , j ^ 
'*’ “l . j , j+2“^^i+l , j+l‘^'^^i+1 ,J ^ 

■ 2^^“i,j+2"^“l,j+l‘*‘^“l.j^^^l+l,j"^i-l,j^ 
■*■ ®1, J ^ ^^1+1 , J+2"*^i+l , J+l'^^^1+1 , j ^ 

“^^l-l,j+2“^^i~l,j+l'*’^^i-l,j^^^ 


(G.2) 


Extrapolations ; 


^1 " ^^i+r^i+2 " ^^i-r^l-2 

^1 " 2^^i+l'*’^l-l^ 

1 (G.3) 

" ■ ■2^^1+2'’^^l+r^^i-l‘*’^l-2^ 

^1 “ “ 6^^i+2"^^i+r^^i-l'’'^l-2^ 

Although all of these expressions, except the last one, are In 
the form which Is traditionally termed "second order accurate", since 
the coordinate transformation Is always chosen so that the 

assignment of an order of accuracy to these expansions becomes ambig- 
uous at best. For example, the central difference first oerlvative ap- 
proximation has as Its leading term in the truncation error the quan- 
tity; - * but, by Itself, this does not give a meaningful ap- 

proximation of the truncation error. There have been some attempts to 
develop estimates of the truncation error in curvilinear coordinates 
which are analogous to the well-known Cartesian error estimates [b,h], 
but more work is needed in this area. 
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APPENDIX H 


Block-Trldlagonal Inversion AlRorlthm 


This algorlthio is basically id -ni ic.il to the one presented by 
Stager [b] in his Appendix 1. Since method is widely available in 

the literature, it is sketched out below purely for the convenience of 
the reader. 

Given the block trldlagonal system of equations: Bu ■■ f, or 

explicitly: 




'*1 


■» M 

^2 ®2 ^2 


“2 


^2 

A3 B3 C3 


“3 

m 

^3 

• • • 


• 


* 

• « • 


• 


• 

• • • 


• 


» 



Vi 


S -1 

1 

1 




/n _ 


(H.l) 


where each A^, and symbol represents a k x k matrix and the vec- 
tors u^ and f j are k-dlmenslonal, the solution may be obtained by a sys- 
tematic application of Gaussian elimination without pivoting. The 
solution scheme can most clearly be developed by splitting the matrix B 
up Into a product of a lower block triangular matrix L and an upper block 
triangular matrix U of the following form: 
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Aj Gj 


A3 G3 


I U. 


I U, 




I IL 


and then finding the solution to the set of matrix equations: 
Lr - f 
Uu - r 


This can be efficiently accomplished by use of the following recursive 
algorithm: 

Forward Sweep 

1) U, - Bj'cj 
r, - Bjlf, 

2) for 1 - 2 to 1 ■ N-1 

‘^l - »i ■ Vl-1 

(H.A) 

'1 ■ - Vi-1> 


Back Sweep 

S. ■ “n - Vn-1 

“n ■ °N ■ 'Wl* 
2) for 1 ■ N-1 to 1 = 1 


“i ■ '1 - "i“i+i 
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As Stagsr points out in ths prsvlously cltsd rsfsrsncs, this 
slgorltha can scconoiodsts ths addition of ssvsral stray matrlcas off of 
ths trldlagonsl strip of ths matrix B, but only at ths cost of mors 
than doubling ths opsratlon count. Also, following ths rsconmsndation 
of Stagsr, ths k x k (hsrs 4x4) matrix invsrsion involved in the above 
solution achsms can bs sfflclsntly calculated by application of the 
same systemisad Gaussian sllminatlon scheme as used in this algorithm. 

Suppose the matrix equation GS M is to be solved for S, where 
G, S, and M are all 4x4 matrices. Let A and B be lower and upper tri- 
angular matrices, respectively, so that 


*11 

0 

0 

0 “ 


T 

^>12 

^3 


®21 

•22 

0 

0 


0 

1 

^23 

"24 

*31 

*32 

*33 

0 


0 

0 

1 

‘34 

_*41 

*42 

*43 

*44 


0 

0 

0 

1 


Then, similar to equation (H.3), the set of equations 
AR - M 

(H.6) 

BS - R 

can be solved by the following algorithm: 

®11 “ ®11 
®21 " 821 
®31 “ *31 
®41 • *41 
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*22 " *22 ’ '* 12*21 
*32 ■ *32 " **12*31 
*42 ■ *42 “ **12*41 

**13 " *13'^*11 

**23 " ^*23 “ **13*21^^*22 

*33 “ *33 “ **13*31 ■ **23®32 
*43 " *43 ' **13*41 " **23*42 

**14 " *14/*11 

**24 ■ ^*24 ■ **14*2p^*22 

**34 " ^*34 ” **14*31 " **24‘"32^^*33 

®44 “ *44 " **14*41 " **24*42 " **34*43 
Repeat for each column of M matrix (j«ltoj"41n this case ) ; 
*^lj “ "ij/*!! 

*"2j " ^"*2j " *21*'lj^/*22 

rsj - (m3j - a3jrjj - 

*^4j “ ^"‘4j ■ *41*^1 j ■ *42^2j ■ *43*^3j^^®44 
*4j “ *^4j 

*3j ‘ '3j - **34*4j 

®2j “ '^*2j " **24*4j " **23®3j 



APPENDIX I 


Body Force Modlflcntloni of Algorithm 


Incluclon of • body force term Is readily done by adding the vector 

0 


P8 


1 




(I.l) 


pugj + pvg2 

to the right-hand-side of equation (3.1) and the corresponding vector 
G* • G/J to the right-hand-side of equation (3.15). Where and g 2 
are nondimenslonal acceleration parameters whose effect is to acceler- 
ate the flow field in the positive x- and y-coordinate directions, or 
equivalently, to accelerate the body and its attached coordinates In 
the negative x- and y-coordinate directions. The generalized time- 
differenced Navier-Stokes equations (4.4) are then accordingly modi- 
fied by adding the term 




( 1 . 2 ) 


where the superscript has been dropped. The Jacobian matrix cor- 
responding to the body force vector G is: 


M 


8q 


0 0 0 0 
gj 0 0 0 

g2 0 0 0 

0 gj g2 ° 


( 1 . 3 ) 


So if, is dlscuacod in Appendix E in regard to the coordinate deriv- 
atives, the acceleration parameters and g 2 are evaluated at time 
(nf6)At, then equation (4.8) is modified by adding G*' to the right 
hand aide and - M^®Aq" to the left hand side. The final algorithm 
expressed in equations (4.9) can now be adjusted to take the following 
form: 

(I + i ,n . X ^ . ^gVrhS(4.8) 

(I. 4a) 

^ '^2 • (l-M"”^®)Aq" (1.4b) 

_n+l n , . n t \ 

q - q + Aq (I. 4c) 

where, it should be emphasized that the acceleration factors and g 2 
are to be evaluated at time (n4d)At In the vector G^. 

The last set of equations may not seem to directly follow from the 
modified form of equation (4.8); that It does, and to the same order of 
accuracy, can be seen from the following discussion. Consider the 
symbolic form of equation (4.8) as modified by the body force term: 

[l + At(M+F+Q)]Aq - R (1.5) 


where M, P, Q, and R represent the body force, ^-derivative , n-deriv- 
atlve, and right-hand-side terms, respectively. Now, It was considered 
to be desirable to factor equation (1.5) In such a way that the body 
force matrix M would have an equivalent effect on the Implicit matrices 
for both the C~ and n-directions. This can be obtained if equation 
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(1.5) Is multlpllsd by ths mstrix S > (I <f dtM) , rssultlng In ths 
squstlont 

[l + At(SP + SQ)]Aq ■ SR (1.6) 

which CSC' thsn bs fsctorsd Into ths sst of squstlons 

(I + AtSP)A^ - SR (I. 7s) 


(I + AtSQ)Aq » TSq (1. 7b) 

to ordsr of At^. Finslly, multiplying these equations by S~^ results 
in: 


[l + At(M + P)]^ - R (1. 8a) 

[l + At(M + Q)]Aq - (I + AtM)M (1. 8b) 

But these last equations are seen to be the symbolic form of equation 
(I. 4a and b). 

As a final note, although the author was unable to prove this, it 
was felt that the magnitude of the acceleration parameters should be 
such that strong diagonal dominance of the implicit matrices is pre- 
served. This results in the criteria 



which tended to be confirmed by numerical experimentation. Note that 
this criteria, if correctly conceived, means that large accelerations 
require small time steps and goes some way toward explaining the great 
difficulties encountered in impulsive starts. 
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APPENDIX J 


Body Boundary Valu« Fonnuiation 


In dotorminlng feho natrlcoi of oquatlon (4.13) it !• convonlcnt to 
consldor eho voctor of dopondont vorinbloo, 4q^ ^ , a componant at a 
tlma. Tha flrat componant, 4p^ ^ , was simply sxtrapolatsd along llnss 
of C to obtain tha "sacond ordar accurato” axprasslon: 


‘“J.! ■ 

(Nota that p and a In this Appandix should proparly ba written as p/J 
and a/J). This boundary condition on density was felt to be an accurate 
approximation since the n-llne spacing was always very rmall near the 
solid bodies studied in this research. The x- and y-momentum differen- 
ces can be exactly re-expressed as: 


./ \0 n+1 . n . n , n 

4(pu)j_, - Uj , 4Pj + Pj j 

A/ \n n+1 , n . n , n 

''‘‘”'>1.1 ■ ''i.i ''‘’i.i ‘’i.i ‘''i.i 


Similarly, the total energy can be exactly expressed as: 


. n /Cvn+l . n . n 

^®1,1 " ^p^i,l ^^i.l ^1,1 


(J.2) 

(J.3) 


(J.4) 


where for the case in which temperature on the body surface is speci- 
fied, the quantities e/p can be evaluated from the equation of state 
(3.5). For all other thermal boundary conditions, equation (J.4) can be 
re-expressed as: 



n 


1 1 + Pi 1 ^Oi'i ■*' 

PX|I X«I ^9^ 


(J.5) 
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Now on aubitltutlng roUclon (J.l) into tha othar aquatlona, tha 
following fora for tha H and K matrlcaa and tha L vactor of aquation 
(4.13) davalopat 


H - 



-10 0 
-u 0 0 
-V 0 0 
-E 0 0 




whom 

E - 1 . It + i(u2 + v2) (J.7) 

and E and AE ate lagged In time unless the surface temperature has been 
sraclfled. 

Although, in a sense, this treatment of the implicit boundary con- 
ditions is rather cavalier, it is this author's experience that the 
stability and, in particular, the accuracy of the algorithm is not very 
sensitive to these boundary conditions (that's not to say that the 
stability and accuracy can not be affected by inappropriate choice of 
the boundary condition). In fact, Steger [b,g] has had good success 
with merely lagging the entire difference vector, Aq , in time. 
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